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Complex systems are often modeled as Boolean networks in attempts to capture their logical 
structure and reveal its dynamical consequences. Approximating the dynamics of continuous vari- 
ables by discrete values and Boolean logic gates may, however, introduce dynamical possibilities that 
are not accessible to the original system. We show that large random networks of variables coupled 
through continuous transfer functions often fail to exhibit the complex dynamics of corresponding 
Boolean models in the disordered (chaotic) regime, even when each individual function appears to 
be a good candidate for Boolean idealization. A suitably modified Boolean theory explains the be- 
havior of systems in which information does not propagate faithfully down certain chains of nodes. 
Model networks incorporating calculated or directly measured transfer functions reported in the 
Uterature on transcriptional regulation of genes are described by the modified theory. 
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I. INTRODUCTION 

Natural systems often involve many types of elements 
interacting in a complicated fashion. The interactions 
may be difficult to describe, and may be mediated in 
ways that are poorly understood. In this situation, it 
is necessary to find a model that captures most of the 
salient features of the system without attempting to de- 
scribe all the details. 

Boolean networks are often constructed to model the 
logic of systems with a complex set of interactions 
[Hi S In this ideaUzation, continuous variables are 
modeled by binary states, and interactions are modeled 
by Boolean update rules. A binary representation is a 
natural approximation for systems whose elements tend 
to take distinct high and low values with sharp transi- 
tions between states. However, even when individual el- 
ements are good candidates for Boolean modeling, qual- 
itative discrepancies between the dynamics of the under- 
lying system and its Boolean idcahzation can arise. 

The relation between continuous and Boolean systems 
has been the subject of some study. In previous work, 
we identified the features that cause a discrepancy in 
the attractor dynamics of small systems [1]. Davidich 
and Bornholdt have shown that through a careful exam- 
ination of the attractor dynamics, a Boolean model can 
be constructed that faithfully reproduces the temporal 
sequence of states obtained from direct binarization of 
a given continuous system. Glass et. al. have ex- 

tensively studied a class of large networks governed by 
piecewise linear differential equations, which involve both 
Boolean and continuous variables. They have noted that 
artifacts introduced by synchronous update contribute 
significantly to the size of the attractor set, and that 
both periodic and chaotic dynamics can be observed, al- 
though chaos is quite uncommon in random networks 
with a connectivity of A: = 2 [1, 0] ■ Magnasco has shown 
that it is possible to construct continuous systems that 
implement any specified Boolean computation Q. The 



present work addresses new dynamical features that arise 
in large networks of generic elements with sigmoidal re- 
sponse functions and explores the extent to which they 
can be understood using Boolean models. 

We consider an illustrative class of continuous mod- 
els and show that information propagation along chains 
plays a key role in determining the qualitative dynamical 
behavior of large random networks. We develop a mod- 
ified Boolean theory incorporating the effects of signal 
decay on certain chains that explains key features of the 
continuous dynamics. Applying the theory to the well- 
known transition between order and disorder in random 
networks [l^ reveals that signal decay has little ef- 
fect on ordered dyanmics but can lead to a substantial 
suppression of disorder. Finally, we study cis-rcgulatory 
functions from the quantitative biology literature and 
show that our theory accounts well for the dynamics of 
random networks constructed using those functions. The 
differences in collective behavior of the continuous sys- 
tems and naive Boolean models are not simple extensions 
of the differences in attractor structure noted in earlier 
studies [3|- They suggest both that caution should be 
taken when making inferences based on Boolean model- 
ing of individual nodes and that appropriately modified 
Boolean models can still provide useful insights. 



II. A CLASS OF CONTINUOUS MODELS 

We study continuous systems with variables {xi} and 
time evolution equations 

Xi{t) = fi[xj{t - nj), Xk{t - r^k)] -^iXi{t) , (1) 



fi{xj,xk) = m 



1 + djX^j + dkx\-\- djkXjXk 
1 + bjx'^j + bkxl + b. 



'jkXjXk 



(2) 



where rji, 7^, the 6's and the o?'s are constant coefficients 
and Tij is a constant time delay associated with the con- 
version of the output of node j into its active form and/or 
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FIG. 1: The transfer functions gi{x,y) (left) and g2(x,y) 
(right). Both functions have the same Boolean idealization. 

a propagation time between node j and its target node i. 
Each node i receives exactly two inputs from randomly 
selected nodes in the network and responds as determined 
by its transfer function fi. The form of the differential 
equations is motivated by studies of genetic regulatory 
networks, in which case the variables represent mRNA 
concentration [Tl| . 

Time delays are included for two complementary rea- 
sons. First, the interactions between elements in the 
systems of interest typically involve a series of events 
that take time to complete. For transcriptional networks, 
the activation or repression of a target gene due to the 
buildup of a particular mRNA requires the translation of 
the mRNA, the folding of the protein, and the transport 
of the protein to the nucleus. In this case, the delays in 
our model equations represent the time required for pro- 
tein levels to build up, along with whatever other post- 
translational processes are required and may be thought 
of as capturing the dominant effect of a set of explicit 
equations for additional variables. Second, the delays are 
necessary to produce simple oscillatory behavior in small 
negative feedback loops. A self-rcpressor in a Boolean 
system produces oscillations rather than fixed points, but 
a continuous self-repressor will not oscillate in the ab- 
sence of a time delay. To get a meaningful comparison of 
large Boolean network models and the underlying contin- 
uous systems, we want to study continuous models capa- 
ble of exhibiting the oscillatory behaviors of very simple 
negative feedback loops. Though there may be cases in 
which the function of negative feedback in a real system 
is just to regulate the level of a fixed point, such feed- 
back would be irrelevant from the perspective of Boolean 
modeling. 

A system that nicely illustrates the role of time delays 
is the repressilator (a loop of three repressors) studied by 



TABLE I: Parameter values for transfer functions. 
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Elowitz and Leibler [T2]. To model the observed oscilla- 
tions, it was necessary to include separate equations for 
the net rate of production of RNA and proteins. Coupled 
equations for RNA concentrations of the form of Eq. [1] 
with r = yield only a stable fixed point With 
sufficiently large time delays, however, the fixed point 
becomes unstable to an oscillatory attractor; it is not 
necessary to include explicit equations for protein pro- 
duction. 



III. ANALYSIS OF TWO SPECIFIC SYSTEMS 

We consider two random network systems. Si and 5*2, 
that have the same Boolean idealization but turn out to 
behave quite differently. Each node in Sa is assigned 
fi = 9a with probability p and fi ~ h otherwise, where 
ga and h are defined by the parameters listed in Table [H 
gi and g2 are plotted in Fig. [T] Both would be approx- 
imated by a Boolean NIF function, which returns a 1 if 
and only if the first input is 1 and the second input is 
0. h corresponds to a Boolean OR function. The choice 
of NIF and OR functions allows us to tunc through the 
order-disorder transition by varying p. The qualitative 
results are not specific to this choice. 

The continuous dynamics are simulated using the 
fourth order Runge-Kutta method described in [3]. The 
time delays {t^} are set to 1, the decay constants {7^} 
are chosen at random from the interval [0.8,1.2], and the 
normalization constants {rji} are chosen such that gi, g2 
and h all have the same saturation value of 100. The 
results presented here are not sensitive to the values of 
{Tij} or {ji}- We studied networks of size N = 1000. 
For each given distribution of transfer functions, we sim- 
ulated between 15 and 30 networks, with 15 to 30 random 
initial conditions each, each network requiring about one 
hour of computation time on a desktop computer. 

To compare the dynamics of the continuous systems 
with their Boolean counterparts, we binarizc the continu- 
ous time series, setting a node's value to 1 (0) if it is above 
(below) a specified threshold. For a given network, let cpi 
be the mean value of node i over time, obtained after 
concatenating time series of equal duration from attrac- 
tors generated by different initial conditions. We focus 
on two quantities: ^ = Avg[ipi], and = Avg [ipi — iff\ , 
where Avg [. . .] denotes an average over nodes, cr^ is the 
average of the variances for a system with binary values 
and 1. Note that frozen nodes, which go to the same 
static value at long times on all attractors, contribute 
zero to (7^ . The concatenation of attractors allows nodes 
that are constant on any single attractor. but not at the 
same value on all attractors, to contribute to . 

Fig. [5] shows the ensemble averages {y) and (cr^) as 
functions of p. Si behaves very much like its Boolean 
idealization, represented by the solid curvesjjSg does not. 
We use established techniques [l3, [H, [13, EHl to deter- 
mine that the Boolean system is ordered for p G [0,0.5), 
disordered for p € (0.5, 1) and critical for p = 0.5 and 
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FIG. 2; (a) The mean fraction of ON nodes. The + and A 
symbols correspond to Si and 52, respectively. The solid line 
/i = 1 — p is the prediction for an infinite random Boolean net- 
work with a fraction p of nif gates and 1 — p of OR gates. The 
dashed curves show the theoretical results discussed in the 
text, (b) The total variance corresponding to the systems in 
(a). The solid curve is the average variance from simulations 
of the Boolean model. 



p = 1. In the large system limit, we expect (cr^) = 
in the ordered regime and a continuous transition to 
(ct^) > in the disordered regime. The nonzero value 
of ((7^) at the critical points is a finite size effect. The 
variance of Si is strongly suppressed in the region where 
the Boolean model is disordered. We note that (/x) = 
implies (c'^) = 0, so for this particular system, we can 
explain the suppression of (cr^) by deriving a theory of 
(/i) . The dashed curves explaining the behavior of Si are 
based on the theory discussed below. 

As discussed by Magnasco in Q , a system whose trans- 
fer functions are globally compatible is capable of execut- 
ing logical operations, where global compatibility means 
that all transfer functions have the same two fixed point 
values; i.e., that when every input is held steady at ei- 
ther the high or low value, every output will also take 
one of those two values. We show here a system (S\) 
that behaves statistically like its Boolean analogue, al- 
though the transfer functions for the two node types are 
not tuned to be globally compatible and the attractors 
are not always steady states. Fig. |3] shows a typical time 
series for a continuous node in 5*1 and a Boolean node 
in networks with p = 0.8. Both exhibit the complicated 
behavior associated with the disordered regime. 

Certain features of continuous systems are known to 
strongly affect the number and nature of attractors in 
small systems and simple rings in ways that are not cap- 
tured by Boolean models It is not clear, however, 
whether these features lead to important effects in large 
systems. The close match for (/x) and (cr^) of and its 
Boolean model (Fig. [2]) suggests that they do not, but 
care must be taken to interpret these results. Much of 
the agreement can be attributed to the fact that (/x) and 



FIG. 3: Typical time series for an unfrozen node in the disor- 
dered regime (p = 0.8) in a continuous system, 5*1, and in a 
synchronously updated Boolean system. Both systems have 
attractors with complex dynamics. 



(cr^) are determined primarily by the fraction of frozen 
nodes in the network, a quantity that does not depend 
on the timing or sequence of updating the nodes @ . 

In the ordered regime, the nearly vanishing number of 
nodes that are not frozen leads to very low values of (c^) . 
The agreement in the disordered regime is less trivial. 
Consider, for example, a ring of nodes containing a ran- 
dom mixture of copiers and inverters. For synchronously 
updated Boolean dynamics, every attractor of a given 
network has a partner in which the values and 1 are 
exchanged for all time steps, and these two attractors 
have basins of the same size. (When the number of in- 
verters is odd, every attractor is its own partner.) Thus 
/i is always 1/2 and is always the maximum value of 
1/4. As discussed in the breaking of on-off symmetry 
in generic continuous systems leads to a collapse of al- 
most all attractors, leaving only two fixed point states for 
rings with an even number of inverters and one oscillat- 
ing state for rings with an odd number of inverters. The 
odd case produces (/x) « 0.5 and (cr^) w 0.25. These are 
not strict equalities because, unlike for the synchronous 
Boolean ring, a node in the corresponding continuous 
system need not spend precisely the same time in the 
high and low states, but the symmetry breaking effect 
is typically small. 0]. In the even case, on the other 
hand, the two fixed point states can have dramatically 
different size basins and therefore exhibit a substantially 
reduced value of . For a ring of size iV = 100 with a 
random mixture of continuous copiers and inverters, we 
find (/x) sa 0.5 and (ct^) ~ 0.12. So in the case of simple 
rings, where (^) and (c^) are determined by dynami- 
cal properties of the attractors rather than numbers of 
frozen nodes, we do not see agreement similar to that of 
S\ and the Boolean system shown in Fig. [51 It appears 
that the complex network of connections between active 
rings in the disordered regime restores the agreement. 
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FIG. 4: (a) gi{x,e) and (b) g2(x,e), with e= 1.12. 



Fig. [3] suggests that Boolean and continuous attractors 
in the disordered regime have similar temporal features, 
but a full characterization of the dynamics of individual 
attractors is beyond our present scope. 

We now turn to the analysis of 5*2, a case where the 
continuous dynamics show substantial deviations from 
the naive Boolean expectations. To understand why 
Si and 5*2 behave differently, we examine the functions 
gi{x,y) and g2{x,y) when the second input is held at a 
low value e = 1.12, the low stable fixed point of h{x, x). 
A Boolean nif function whose second input is held at 

acts as a copier on its first input, so we expect the 
functions gi{x^e) and g2{x,e) to behave like copiers. As 
shown in Fig. 21 gi{x, e) has two stable fixed points, but 
g2{x, e) has only one stable fixed point. 

As noted in , the faithful propagation of information 
along chains of nodes requires two stable fixed points in 
the transfer function. We refer to the loss of information 
along chains of nodes lacking a second stable fixed point 
as "propagation failure." Because g2 does not have the 
required fixed point structure, propagation failure may 
cause the system S2 to behave differently from its naive 
Boolean idealization. We now present a modified Boolean 
model that accounts for propagation failure and agrees 
well with simulations, as shown by the dashed curves 
in Fig. [2] The success of this theory indicates that for 
large random networks, propagation failure is the pri- 
mary source of the measured discrepancy between the 
continuous system and the original Boolean model. 

A modified Boolean theory for S2 is derived by noting 
that some of the nominal nif nodes are actually better 
modeled by the Boolean off function. Consider a chain 
of continuous nodes where each has the transfer function 
g2{x, y), and where node i is the first input (x) into node 

1 + 1. Because g2 lacks the necessary high fixed point, 
propagation failure prevents nodes far down the chain 
from ever rising above threshold in response to a high 
input signal to the first node in the chain (i = 1). Though 
they may be initialized at a high value, they will always 



stay low after a brief transient. Let node m be the first 
node in the chain that cannot rise above threshold due 
to propagation failure. We will model node i with a nif 
function if i < to, and with an OFF function if i > to. 
The value of to for the a particular chain of 52 nodes will 
depend on the parameters in g2, the choice of threshold, 
and the value of the high input into the first node of the 
chain. Approximating the system as a random Boolean 
network with fractions r, g, and 1 — r — (/ of OFF, nif, 
and OR nodes, respectively (neglecting correlations in the 
positions of the OFF nodes), (/i) can be calculated exactly 
as the stable fixed point of the bias map [l^ 

Pt+i =q{pt' Pt) + {'^-r- q) {2pt ~ pi) , (3) 

where pt represents the average fraction of nodes with a 
value of 1 at time t. The stable fixed point is 

(Ai) =max|o,^^^| . (4) 

Let n be the value of to for iS'2, computed using the 
average value of the high input signals that arise from 
the dynamics. For a large random network, the fraction 
of nodes that are assigned the OFF function is r = . so 
the fraction of nodes that truly act as nif \s q— p — p"' . 
Substitution into Eq. [4] gives 

(,.)=max|o,i^^}. (5) 

The prescription for finding off nodes does not explicitly 
account for the failure of propagation around rings of g2 
nodes smaller than n, but such loops are very rare in 
random networks. 

The dashed curves in Fig. [5] show {p) for n G {2,3,4}, 
with n increasing to the right . The simulations suggest 
a crossover from n — A for p < 0.5 to n = 3 for p > 
0.5. The reduction of n arises because the average value 
taken by nodes that are above threshold decreases as p 
increases, which implies lower input values to chains of g2 
nodes. The switch from ri = 4 to n = 3 associated with 
the function 32(2;, e) of Fig. 2] and our chosen threshold 
of 10 occurs when the initial input to a chain is about 
58. Simulations reveal that the average value of nodes 
above threshold crosses 58 at p w 0.51, which corresponds 
reasonably well with the crossover observable in Fig. O 

The suppression of disorder, as indicated by small val- 
ues of (cr^) for p > 0.5 in Fig.[2)^b), is caused by the effec- 
tive insensitivity of nodes at the end of sufficiently long 
chains. Because nodes in this set have approximately 
fixed values, all nodes receiving both inputs from this set 
will also have approximately fixed values, leading to a 
cascade of effectively frozen nodes. The net result is a 
substantial loss of sensitivity for the network as a whole. 

Further insight into the order-disorder transition 
comes from examining the trajectory traced in the q — r 
plane as p increases from to 1, shown in Fig. O The 
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FIG. 5: (a) Trajectory in the q — r plane, as predicted 
by the mean field model, for fixed decay lengths n G 
{2,3,4,5,10,20}. The curves go from low to high n when 
viewed left to right, (b) Fit of actual data with varying n. 
See text for details. 



shaded sector is the disordered regime of the Boolean sys- 
tem, where the slope of the Derrida plot exceeds unity 
(lOl . [iGj . The unshaded sectors correspond to ordered 
regimes, and the boundaries to critical systems. Fig.[5l^a) 
shows trajectories corresponding to different n. We can 
measure (/z) in simulations of ^2 and determine r and 
q from the relations r -\- q = p and Eq. 3) The results 
are shown in Fig. [Hb). Only the results for (/z) > are 
plotted; the theory cannot determine unique values for r 
and q in the upper unshaded sector because every (g, r) 
pair gives (/z) = 0. This plot again reveals a shift from 
n = 4 to n = 3, but also shows that the system skirts 
a critical boundary as p is varied. The dynamical sup- 
pression of disorder in this case provides a mechanism 
for keeping a network in the critical regime over a wide 
range of parameter values. We note, however, that this 
is not necessarily a generic effect. Functions producing 
larger values of n would permit some degree of disorder. 



IV. RESULTS FOR REPORTED TRANSFER 
FUNCTIONS 

It is instructive to examine suggested or measured 
transfer functions from the literature to see whether col- 
lections of similar functions would faithfully execute their 
nominal Boolean logic. In Boltzmann weights are 
used to compute the probability that RNA-polymerase 
will bind to the promoter region of a gene as a function 
of the concentrations of the transcription factors for that 
gene. If we take the transfer functions to be proportional 
to the probability functions, we find that they lack the 
necessary fixed points for propagating information and 
that random networks built from them exhibit strong 
suppression of disorder. We note, however, that there is 
evidence for cooperative effects and post-transcriptional 
processes that influence the effective transfer function, 
so we do not necessarily expect the systems considered 
below to be representative of real regulatory networks. 

Consider a continuous system with two transfer func- 
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FIG. 6: Mean node value of a network using nand and OR 
transfer functions from and a threshold that is (-1- sym- 
bols) below, and (A symbols) above the intermediate fixed 
point defined in the text, p is the fraction of nand nodes. 
Our mean field theory (dashed curves) gives a much better 
approximation of the dynamics than does the naive Boolean 
model (solid curve). 



tions defined in as implementing nand and OR logic. 
A mean field theory analogous to Eq. O can be con- 
structed. The results in this case depend upon the choice 
of threshold for the binarization of the time series. Two 
plots for (/i), using different threshold values, are given 
in Fig. [SI The solid curve shows the behavior of the 
straightforward Boolean idealizations, and the dashed 
curves show the predictions made by our mean field the- 
ory. The theory predicts ^cr^) ~ for all fractions p 
of NAND nodes because propagation failure pushes the 
system into the ordered phase, an effect confirmed by 
simulations. As above, accounting for the inscnsitivity 
of some nodes allows a reasonably accurate prediction of 
the dynamics of the continuous system. 

In a random Boolean network of nand and OR gates, 
the fixed point in the bias map becomes unstable to a 
2-cycle at p = V3/2. Oscillations in bias occur also in 
continuous systems with long time delays and short decay 
times, but not if these time scales are comparable. In 
the latter case, the nand nodes find stable intermediate 
values. Fig.[6]shows the dynamics with the choice = 1 
and 7i G [0.8, 1.2] with two different threshold values, one 
above and one below the intermediate fixed point of the 
all NAND system. Analysis of continuous systems with 
oscillating bias is beyond the scope of this work. 

Ref. provides no obvious guide for choosing the 
normalization constants for the transfer functions. For 
Fig. [5] wc used rji = ri2 — 1000. We tried every combi- 
nation of (?7i, 7^2) with m e {200, 400, 600, 800, 1000} and 
found no pair that produced behavior significantly differ- 
ent from those shown. In this case, as in most cases, the 
NAND nodes took clearly separated high and low values 
for most values of p, but the OR nodes typically converged 
to a single value due to the absence of a low fixed point. 

A second example of a biological transfer function is 
the lac cis-regulatory input function studied in [l8| . The 
lac promoter requires both IPTG and cAMP to be pro- 
duced. Setty et. al. discuss both the real input function 
and an idealized continuous and function. Neither of 
these functions satisfies the criterion for faithful infor- 
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mation propagation in a network. For high values of the 
first input, the function does not have a stable high fixed 
point for the second input. A large network of such ele- 
ments would not express its nominal Boolean logic. 

It has been shown both theoretically and experimen- 
tally that si gna ling cascades can create switch-like re- 
sponses [l^, |20| . w\ in genes several steps downstream 
from the initial stimulus. In transcriptional regulation, 
one mechanism for achieving the required characteristics 
of the transfer functions involved was proposed in [2^ . 
wherein an evolutionary model generated steeper trans- 
fer functions by introducing auxiliary binding sites for a 
transcription factor. Another mechanism has been pro- 
posed in recent work by Buchler et. al., who emphasize 
the possible role of protein or RNA sequestration in cre- 
ating sharp transfer functions [1^, [1^ . We note here that 
the demonstration of an effective Hill function with high 
cooperativity does not by itself imply that a cascade of 
similar elements would successfully propagate informa- 
tion. One needs to know how to normalize the output 
levels to see whether there is a high fixed point, and this 
requires knowledge of the level needed to activate the 
downstream target. Nevertheless, the use of cooperativ- 
ity and sequestration to construct elements with nearly 
Boolean behavior appears promising. Our present work 
highlights another relevant issue for synthetic or natural 
transcriptional circuits: the interactions of transcription 
factors may make it difficult to simultaneously meet all 
of the requirements for signal propagation for genes with 
two or more inputs. Attempts to experimentally demon- 
strate a sharp, two-input logic gate may yield new in- 
sights into the feasibility of implementing complex logic 
in biological systems. 

It is not clear, however, that biological networks need 
to be able to carry out arbitrarily complicated logical 
operations. Our results show that faithful propagation 
down long chains is not crucial for implementing the 
Boolean dynamics arising in the ordered and critical 
regimes in systems of a few thousand elements. More- 



over, if operation near criticality is advantageous for a 
biological system, as suggested in [l| , suppression of dis- 
order may be beneficial. We have identified a dynamical 
mechanism for accomplishing this, which may be a useful 
alternative to alteration of the network architecture. A 
limited range of propagation could actually be an impor- 
tant feature of real biological networks. 



V. CONCLUSION 

In conclusion, we have seen that the detailed form of 
continuous transfer functions can have a qualitative ef- 
fect on the dynamics of large random networks. If all 
of the transfer functions have suitable fixed point struc- 
tures such that signals can propagate down chains of ar- 
bitrary length, fundamental statistics of the attractors 
arc similar. This indicates that disordered systems are 
not prone to the type of attractor collapse that occurs 
on simple rings. When some of the transfer functions do 
not have a suitable fixed point structure, the dynamics in 
the disordered regime are strongly affected. In the case 
studied above, the suppression of disorder leads to an ex- 
tended domain in the parameter space where the system 
is close to a critical Boolean network. A successful the- 
ory of these effects has been obtained within a Boolean 
framework by accounting for the failure of information to 
propagate down long chains of nodes. Further character- 
ization of the temporal structure of individual attractors 
in the disordered regime is needed to fully understand 
the importance of these phenomena in large networks. 
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